// An OpenCL implementation of a multithreaded Mersenne Twister PRNG.
 
#define MT_MM     9
#define MT_NN     19
#define MT_WMASK  0xFFFFFFFFU
#define MT_UMASK  0xFFFFFFFEU
#define MT_LMASK  0x00000001U
#define MT_SHIFT0 12
#define MT_SHIFTB 7
#define MT_SHIFTC 15
#define MT_SHIFT1 18

typedef struct {
  unsigned int matrix_a;
  unsigned int mask_b;
  unsigned int mask_c;
  unsigned int seed;
} MersenneTwisterConfig;

typedef struct {
  unsigned int mt[MT_NN];
  int iState;
  unsigned int mti1;
} MersenneTwisterState;

// Initialises the state structure and must be called before any calls to
// MersenneTwisterGenerate*() are made.
void MersenneTwisterInitialise(MersenneTwisterState *const state, __local const MersenneTwisterConfig *const initialConfig)
{
  state->mt[0] = initialConfig->seed;
  for(int i = 1; i < MT_NN; ++ i) 
  {
    state->mt[i] = (1812433253U * (state->mt[i-1] ^ (state->mt[i-1] >> 30)) + i) & MT_WMASK;
  }

  state->iState = 0;
  state->mti1 = state->mt[0];
}

// The main library function which any thread can call any time to receive the
// next integer random number (see below for floating-point variant) from the
// Twister referenced in the state structure.
unsigned int MersenneTwisterGenerateUInt(MersenneTwisterState *const state, __local const MersenneTwisterConfig *const initialConfig)
{
  const size_t tid = get_global_id(0);

	int iState1 = state->iState + 1;
	int iStateM = state->iState + MT_MM;

	if(iState1 >= MT_NN) iState1 -= MT_NN;
	if(iStateM >= MT_NN) iStateM -= MT_NN;

	unsigned int mti = state->mti1;
	state->mti1 = state->mt[iState1];
	unsigned int mtiM = state->mt[iStateM];

	unsigned int x = (mti & MT_UMASK) | (state->mti1 & MT_LMASK);
	x = mtiM ^ (x >> 1) ^ ((x & 1) ? initialConfig->matrix_a : 0);
	state->mt[state->iState] = x;
	state->iState = iState1;

	// Tempering transformation.
	x ^= (x >> MT_SHIFT0);
	x ^= (x << MT_SHIFTB) & initialConfig->mask_b;
	x ^= (x << MT_SHIFTC) & initialConfig->mask_c;
	x ^= (x >> MT_SHIFT1);

	return x;
}

// Floating-point PRNG is a thin wrapper over integer PRNG, which MTs work in.
float MersenneTwisterGenerateFloat(MersenneTwisterState *const state, __local const MersenneTwisterConfig *const initialConfig)
{
  return (float)(MersenneTwisterGenerateUInt(state, initialConfig)) / 4294967295.0f;
}

__kernel void TestMersenneTwisters(__global const MersenneTwisterConfig *const initialConfigs, __global float *const outArr, const int nRandomNumbers)
{
  __local MersenneTwisterConfig sMem[512];
  __local MersenneTwisterConfig *const initialConfig = &sMem[get_local_id(0)];

  // Copy our initial configuration into local memory to reduce read times.
  initialConfig->matrix_a = initialConfigs[get_global_id(0)].matrix_a;
  initialConfig->mask_b   = initialConfigs[get_global_id(0)].mask_b;
  initialConfig->mask_c   = initialConfigs[get_global_id(0)].mask_c;
  initialConfig->seed     = initialConfigs[get_global_id(0)].seed;

  // Our state structure is too large to fit into shared memory and must sit
  // in slower private memory.
  __private MersenneTwisterState mtState;
  MersenneTwisterInitialise(&mtState, initialConfig);

  // "Warm up" the Twister to avoid initial correlation with other Twisters.
  for(int i = 0; i < 10000; ++ i) {
    MersenneTwisterGenerateUInt(&mtState, initialConfig);
  }

  // Generate a subset of the random numbers.
  for(size_t i = get_global_id(0); i < nRandomNumbers; i += get_global_size(0)) {
    outArr[i] = MersenneTwisterGenerateFloat(&mtState, sMem); // &sMem
  }
}

// This kernel is meant to be identical to the above *except* for the PRNG
// generation.
//
// Measuring performance by difference is the most accurate way to gauge the
// the performance of MersenneTwisterGenerateFloat() within the kernel, which
// is the most useful figure to clients of the library functions.
__kernel void TestNoOp(__global const MersenneTwisterConfig *const initialConfigs, __global float *const outArr, const int nRandomNumbers)
{
  __local MersenneTwisterConfig sMem[512];
  __local MersenneTwisterConfig *const initialConfig = &sMem[get_local_id(0)];

  // Copy our initial configuration into local memory to reduce read times.
  initialConfig->matrix_a = initialConfigs[get_global_id(0)].matrix_a;
  initialConfig->mask_b   = initialConfigs[get_global_id(0)].mask_b;
  initialConfig->mask_c   = initialConfigs[get_global_id(0)].mask_c;
  initialConfig->seed     = initialConfigs[get_global_id(0)].seed;

  // Our state structure is too large to fit into shared memory and must sit
  // in slower private memory.
  __private MersenneTwisterState mtState;
  MersenneTwisterInitialise(&mtState, initialConfig);

  // "Warm up" the Twister to avoid initial correlation with other Twisters.
  for(int i = 0; i < 10000; ++ i) 
  {
    MersenneTwisterGenerateUInt(&mtState, initialConfig);
  }

  // Generate dummy writes inbto the random memory buffer.
  for(size_t i = get_global_id(0); i < nRandomNumbers; i += get_global_size(0)) 
  {
    outArr[i] = 0.0f;
  }
}